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We investigate a multilayer stack of dipolar Bose-Einstein condensates in terms of a simple Gaus- 
sian variational ansatz and demonstrate that this arrangement is characterized by the existence 
several stationary states. Using a Hamiltonian picture we show that in an excited stack there is 
a coupled motion of the individual condensates by which they exchange energy. We find that for 
high excitations the interaction between the single condensates can induce the collapse of one of 
them. We furthermore demonstrate that one collapse in the stack can force other collapses, too. 
We discuss the possibility of experimentally observing the coupled motion and the relevance of the 
variational results found to full numerical investigations. 
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I. INTRODUCTION 



Over the past fifteen years Bose-Einstein condensates 
(BECs) have become an active field of theoretical and ex- 
perimental investigations. The experimental realization 
of dipolar BECs with chromium atoms [Ij (for a recent 
review and further references see Ref. [5]) has opened the 
way to actually observe the theoretically predicted phe- 
nomena of radial and angular rotons, anisotropic solitons, 
biconcave shapes of the ground state, etc. |3:,^5 which 
should exist in such condensates with an additional long- 
range interaction. Moreover, dipolar BECs are of special 
importance since the interactions between the atoms can 
be tuned from predominantly short-range to the domi- 
nance of the long-range dipole-dipole interaction (DDI) 
by manipulating the s~wave scattering length via Fesh- 
bach resonances. 



Recently, multilayer stacks of dipolar BECs have been 
in the focus of theoretical investigations jH [7], even 
though they have not yet been realized experimentally. 
In Refs. |6l numerically exact calculations on grids 
have been performed via imaginary time evolution, re- 
vealing structured ground state wave functions and the 
roton instability. By contrast, the purpose of this pa- 
per is to investigate such a multilayer stack of dipolar 
BECs in the framework of a variational approach with a 
Gaussian type orbital for each layer. In this way we shall 
not be able to catch exotic features of dipolar BECs such 
as structured wave functions or the roton instability but 
we can determine different stationary states and study 
in particular the dynamics of the coupled BECs. Fur- 
thermore, we show that excited BECs exchange energy 
and that the scattering length as well as the distance be- 
tween the BECs have a strong influence on this energy 
exchange. For highly excited BECs this energy exchange 
can induce the collapse of one of the BECs, and the col- 
lapse of a single BEC in the stack can cause other col- 
lapses, too. 




FIG. 1. A multilayer stack of dipolar BECs, each placed in a 
trap that is assumed to be very oblate. The single condensates 
are displaced by a distance A, respectively, and coupled by 
the long-range dipole-dipole interaction (cf. [3). 



II. THEORY 

In accordance with Rcfs. [B', '7|, we investigate a stack 
of Ns dipolar BECs. Each of the condensates is arranged 
in an axisymmetric trap that is assumed to be very oblate 
(see Fig. [T]) and the traps (and consequently the BECs) 
are displaced by a distance A from each other in the 
z-direction. 

At very low temperatures the quantum gas in each 
condensate j can be described by a single wave func- 
tion ipj{r,t) whose dynamics obeys the Gross-Pitaevskii 
equation (CPE) 

A + ftrap + Vc + Vd) V^j {r,t)= idt^pj {r,t), 
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where 



Vtrs 



{N, + l- 2j)A 



is the potential energy caused by the harmonic traps each 
condensate is placed in. These traps are arranged in 
the coordinate system in a way that the z — 0-plane 
is always at the center of the stack. Here the GPE is 
written in dimensionless form obtained by introducing 
"atomic" units |B]: Using the mass m of the bosons and 
their magnetic moment /Lt, we define a "dipole length" 
Qd = mnQiJ,^ /{2iTh^), a unit energy Ed — h? /{2ma'\) and 
a unit frequency ujd — Ed/h. The trap geometry is de- 
termined by a mean trap frequency Cj = (WpW^)^/^ and 
a trap aspect ratio A = ujz/ujp. The parameters 7p^z 
in the GPE are connected to the trap frequencies by 



7p,; 



ujp z/{2uJd)- Furthermore, we assume each BEG 



to consist of the same number of particles N so that 
we can apply a particle number scaling r — > Nadr and 
E EdE/N"^ to make the interaction potentials inde- 
pendent of the particle number. The term 

T/e = 87r-|^,(r)|2 
a-d 

represents the contact interaction, which depends on the 
s-wave scattering length a, measured in units of a^, and 



dv — ^ 



r — r' 



1=1 



is the potential energy caused by the dipole-dipole in- 
teraction (DDI) between all the atoms. The long-range 
nature of the DDI demands the summation 1 = 1,... ,Ns 
and couples all the individual condensates. 

We investigate the multilayer stack of dipolar BEGs 
variationally using a cylindrically symmetric Gaussian 
trial wave function 

N, 



where each condensate j is described by a single wave 
function 



.- ^ 1/2 / \ 1/4 




Here ^ = J^^ '^ '^^gj ^"^^ ~ +^^0 ^''^ complex val- 
ued and time-dependent width parameters that are split 
into their real and imaginary parts and have to satisfy 
A^,A,j > 0. The wave functions i^j{r) are assumed to 
be non-overlapping and normalized: 



We have to choose the value of A in such a way that on 
the one hand it is large enough so that we can ensure the 
wave functions to be non-overlapping and on the other 
hand that it is small enough to make the interaction be- 
tween the BEGs sufficiently large. 

To investigate the stack of multilayer BEGs, we will 
first calculate the mean-field energy and then apply a 
time-dependent variational principle. The dynamics of 
the system will be calculated using an equivalent Hamil- 
tonian picture which reveals the physics of the system 
in a more transparent way than that of complex width 
parameters. 



A. Calculation of the mean-field energy 



The mean-field energy E^f of the arrangement is given 

by 

En.i = Jd\ r{r) (^-A + Vt.^p + ^V, + ^Kiip)^(r) 

where, as usual, we have to insert a factor of 1/2 for the 
contact and the dipole-dipole interaction to avoid a dou- 
ble counting of the two-particle interactions. The contri- 
butions of the kinetic energy, the potential energy in the 
harmonic traps and the contact interaction can easily be 
calculated: 



dV V*(r)(-A)^(r) = ^ 2A^ + 2 

j=i 

+ 4r 



. [4/ 



4, 



/dV^*(r)W(r) = 5: — ^+ _ 

\ r A r- 

^ydVV/(r)K^(r-) = E^^4.-V^ 



(1) 
(2) 

(3) 



All integrals occurring here are elementary. The contri- 
bution of the DDI, given by 



1 



dV i/;*(r)T/dV(r), 



is more difficult to calculate. We do this Fourier trans- 
forming the DDI potential 




3„' 



1 - 3 



(2^) 



3/2 ^ -p ) (r-r'Y 



:F{Y^\Mr') 



. 1=1 
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where {...} and J^~^ {...} denote the symmetrical form 
of the Fourier transform and its inverse. In the last step, 
we have applied the convolution theorem. For the first 
Fourier transform one obtains [S] 



(2^) 



3/2 



1-3 



(r- 



47r /3kl 



- 1 



The Fourier transform of J2i {''Pi (''') I can be calculated 



in a straightforward way, and the result is 



T{J2\Mr)\ 




where the term ~i{Ns + 1 — 2l)Akz/2 occurs because of 
the displacement of the individual condensates in the z- 
direction. We can now write the contribution of the DDI 
to the variational mean-field energy as 



-1 J 47r 



3/2 I fc2 



- 1 • exp - 



.{Ns + l- 21) Ak, 



(4) 



where we first integrate over r and afterwards solve the integral over k resulting from the inverse Fourier transform. 
The variational dipole interaction energy finally reads 



dV '^P*{r)Vd^|J{r) 



1 



E 



ds 



■ exp 



4 2 
sin 



■ exp 



45 

2(i - O^A^ 



1 

13/2 



2g5/2 



(5) 



where we have introduced the function g = a^i/S + 

2a/,(s-l) and new variables a//" = [4,-,]" ' + [^^..^ ' 
for brevity. It is useful to split the sum into a part j = I 
and a part j ^ I (i.e. a term that describes the DDI 
between all the atoms belonging to the same condensate 
and a term that describes the interaction between differ- 
ent condensates) because the remaining integral can be 
solved analytically for j = I- In this case the result is 



d'rrir)vjl;'^^ir) 



2 ^^my^j f 2 Syy^ arctan - 1 



(6) 



with rij = Azj/Apj. The mean- field energy of the whole 
stack is then given by subsuming Eqs. ([l} - (6| with the 
constraint j I for the summation in Eq. (|5 1 . 



B. Dynamics of multilayer stacks of dipolar BECs 

To describe the dynamics of the stacked BECs, we ap- 
ply a time-dependent variational principle (TDVP) to 
solve the time-dependent GPE 

with the parameter-dependent wave function ■0(t) — 
ip{A{t)) where A — {ApjAzY is the set of time- 
dependent variational parameters. The TDVP has al- 
ready been applied to condensates with an attractive 
gravity- like 1/r-interaction [HP and to dipolar BECs [S]. 
The application of the McLachlan variational principle 
[TT] leads to a set of first order differential equations, 
describing the time evolution of the real and imaginary 
parts of the width parameters. 

Since it is more descriptive, we present the dynam- 
ics of the arrangement using an equivalent Hamiltonian 
picture, describing the wave function of each condensate 
by a particle moving in a 2A^-dimensional space. The 
Hamiltonian form is obtained introducing generalized co- 
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ordinates 



2^1 



1 



and generalized momenta defined by 



V2A 



Ppj 



With this definition, the generalized coordinates can di- 
rectly be interpreted as the "extension" of the condensate 
in the radial and z-direction, respectively. The Hamilto- 
nian is formally obtained by substituting 

i^a ' ' ' ^3 ) ^ ("^PJ ' ' Pp3 ' P^3 ) 

in the mean-field energy, which results in the Hamiltonian 
H = T + V, where 



T : 



Y.PP3 

3 = 1 



is the kinetic energy and 



rP- An2 ^ Ip'ip3 ^ Iz'dzj^ rz 

, = 1 ^Pi ^^3 V 



1 



E 



ds / 1 



1 2 



3 V TT 



E 



1 



cxp 



la 1 1 1 + 277j — St^j arctan y^rjj — 1 



fld 6 g2 g^^. _ 1 



exp 



45 



(7) 



is the external potential in which the "particle" moves. 
Here the substitution leads to g = (g^; +q'ii)/2 + '^{q^pj + 

The time evolution of the particle in the 2iV5- 
dimensional space {qp,qz) is then determined by the 
Hamiltonian equations of motion 



dH 



and 



OH 



P,Z3 



(8) 



which, after backward substituting {qpj,qzj,Ppj,Pzj) — > 
(A^ , A"^ , , A[,j) , yield the same equations of motion ob- 
tained by applying the TDVP, hence the Hamiltonian 
form is equivalent to describing the condensates using 
parameter-dependent Gaussian trial wave functions. 

Summarizing all the coordinates q = {qp,qzY in one 
vector, the Hamiltonian equations of motion can easily 
be brought into the form 



(9) 



where f{q) — ~2dV/dq is a function of all the coordi- 
nates and the external parameters. Depending on the 
value of a /ad and the number Ng of BECs in the stack, 
V is characterized by one minimum and several saddle 
points, each corresponding to a stationary state of the 
multilayer stack of BECs, given by the fixed points 



0, 



0. 



To investigate the stability of the different states and 
the motion in the very vicinity of their fixed points, we 
follow the usual procedure and linearize Eq. ^ around 
one of its fixed points go , which results in the differential 
equation 



u = Ju. 



(10) 



Here the vector u = q qo denotes the deviation of the 
particle from the fixed point and J = df{q)/dq\^^^^^ is 
the Jacobian matrix. Eq. ( 10 ) is that of a system of 2Ns 



coupled oscillators and can simply be solved using the 
ansatz u = itoe"* with a complex parameter k. Insert- 
ing this ansatz into the differential equation yields the 
eigenvalue equation 



[J -k^)uq = 



(11) 



where are the eigenvalues of the Jacobian matrix. It 
can easily be shown, that J is symmetric, consequently 
the eigenvalues are purely real. 

If all the eigenvalues are negative, u oscillates around 
the fixed point according to it ~ e*'^' with a frequency 
Lu ~ V— and the fixed point is stable. If one of the 
eigenvalues is positive, there is a contribution u ^ e''* 
and the fixed point is unstable. 
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a) -6.30x10" 



FIG. 2. Mean-field energy of the stationary states for a stack 
of = 3 condensates for a distance A = 0.07, scaled trap 
frequency N'^'y = 1300 and aspect ratio A = 340 with Ni, = 6 
different branches. Inset: The different branches arise in tan- 
gent bifurcations at difi'erent values of the scattering length. 
The symbols "s" and "u" denote weather of not the individual 
BECs are stable or unstable. 



III. RESULTS 
A. Stationary states 



s-s-s 
s-u-s 
u-s-s 
u-u-s 
u-s-u 
u-u-u 




b) 5.0x10° 
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We first want to focus on the stationary states of the 
system and show that it is characterized by one stable 
and several unstable stationary states. Even though the 
unstable states cannot be observed experimentally, they 
give a clear picture of the structure of the external po- 
tential V in Eq. ([t]) , which also determines the dynamics 
of the system, discussed in Sec. |IIIB| 

For a given set of physical parameters N'^^p^z, CL/oid, 
A and the stationary states of the system of inter- 
acting dipolar BECs are calculated solving the Hamilto- 
nian equations of motion for vanishing time derivatives 
in Eq. ([s]), i.e. qp,zj = 0, Pp.zj = 0. This results in 
Pp = Pz = ^ and a 27Vj-dimensional, highly nonlinear sys- 
tem of equations for the generalized coordinates qpj , qzj 
which is solved numerically after providing initial values 
for Qpj and qzj (j = 1, . . . , N^). 

For a single EEC there exist two different stationary 
states for a scattering length above the critical value. One 
of them is stable ( "s" ) and one is unstable ( "u" ) [S] . By 
analogy with that we label the states of the multilayer 
stack of BECs by a combination of "s" and "u" mean- 
ing that for A — >■ cxD (i.e. vanishing interaction between 
the single BECs) the stack would be divided into single 
BECs, each in a stable or an unstable stationary state. 

Fig. [2]shows the mean-field energy for stationary states 
of the stacked BECs for Ns — i condensates, a distance 
of A = 0.07 and a trap geometry defined by N"^^ = 1300 
and an aspect ratio of A = 340 (this trap geometry has 
been used for calculations in Ref. [7] and will also be 
used for all calculations in this paper). As can be seen 



FIG. 3. Eigenvalues of the Jacobian matrix J for a stack of 
Ns = ?) BECs and different arrangements in dependence on 
the scattering length. The physical parameters are the same 
used for the calculation of the mean-field energy in Fig. [2] 
The state "s-s-s" (solid red lines) is the only one with ex- 
clusively negative eigenvalues, and in general the eigenvalues 
corresponding to a motion in z-direction (a) differ by several 
orders of magnitude from those corresponding to the motion 
in radial direction (b). 

there are several branches of the mean-field energy. The 
number of the stationary states can be explained assum- 
ing two different stationary states of each condensate ( "s" 
and "u" ) which would result in 2^= possible arrangements 
for the stack. Since some of the arrangements are physi- 
cally equivalent, because they only differ by an inversion 
with respect to the z = 0-plane, the number of inde- 
pendent different arrangements, and hence the number 
of branches, is 




2^.72-1 ^ 2^»-i for N, even 

2W-i)/2^2^=-i for TV, odd. 



The different states arise in tangent bifurcations at dif- 
ferent values of the scattering length. We note that two 
states arising together always belong to the same type of 
symmetry with respect to the z = 0-plane: Either both 
states are symmetric or both are antisymmetric. 

The eigenvalues of the Jacobian matrix for the same 
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set of physical parameters used to calculate the mean- 
field energy are shown in Fig. [3] and the range of the 
scattering length corresponds to that of the inset in Fig. 
[2j Fig. [3^ and Fig. [3]d show the eigenvalues correspond- 
ing to the motion in p- and z-direction, respectively, 
which differ from each other by several orders of magni- 
tude because of the large trap aspect ratio of A = 340. 
The state labeled "s-s-s" is the only one with exclusively 
negative eigenvalues and consequently is the only sta- 
ble state. All the other states have at least one positive 
eigenvalue and are unstable. Note that, for this reason, 
numerically exact calculations [7] can only access the sta- 
ble ground state but not the excited states. By contrast 
all stationary states are independently of their stability 
accessible by the variational approach, which is one of 
the big advantages of that method. 

It is important to note that the stability properties 
of the multilayer stacks do not depend on symmetry as- 
sumptions for the wave function. Since the ansatz with a 
cylindrically symmetrical trial wave functions implies re- 
strictions, we additionally investigated the stability of the 
stationary states using an extended three-dimensional 
trial wave function 



' 93 4' 4' 4» \ 



1/4 



exp I xA^x^ + ■■vAy^y'^ + iA,j ( z ■ ' ^ "^^^^ 



which in principle also allows to study anisotropic trap 
geometries. It turns out that the stability of the station- 
ary states is not affected by this generalization, in par- 
ticular the mode that becomes unstable when one goes 
below the critical scattering length remains cylindrically 
symmetric. Thus the general features found in the two- 
dimensional approach remain valid. 

We also note that the simple variational ansatz con- 
firms the dependence of an increasing critical scattering 
length when one adds more condensates to the stack pre- 
sented in Ref. [7] . Moreover it shows that the wave func- 
tions of the central condensates are accumulated near the 
symmetry axes while the outer BECs are more extended 
in the radial direction. Of course, as mentioned above, 
non-Gaussian structures in the wave function cannot be 
revealed with this approach. 
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FIG. 4. Time evolution of a 2-layer stack of dipolar BECs 
described by a particle moving in the external potential Q 
for A = 0.035 and a scattering length of a/ad = —0.1. At 
t = 0, the particle is placed in the minimum of the potential V 
with nonzero initial momentum ppi ^ 0. The two condensates 
represented by the particle exchange their excitation energy 
periodically. For a small initial momentum ppi (a) the energy 
exchange happens in a shorter period of time than it does for 
large values of the initial momentum (b). 



the Qp-coordinates of the particle) for this situation with 
a small value of the initial momentum ppi . It is calculated 
solving the Hamiltonian equations of motion (|8| using a 
Runge-Kutta algorithm. As can be seen, the extension 
of the excited BEG {j = 1) begins to oscillate around 
its stationary value quickly. Due to the interaction, its 
excitation energy is transferred to the condensate j = 2 
on a larger time scale, and the two BEGs continue ex- 
changing energy periodically. The energy exchange can 
also be observed for more than two BEGs, the difference 
being that the whole energy is not exchanged between 
two single BEGs but transferred to all the others in the 
stack. 

For small excitations the coupled motion of the BEGs 
can be described by linearized equations of motion, which 
reveal normal modes that show a coupling between the 
radial motion of all BEGs and the motion in z-direction, 
respectively, caused by the long-range DDI. An inves- 



B. Dynamics 

We demonstrate the coupled motion of the BEGs by 
placing the particle representing the condensate wave 
functions at the stable fixed point of the Hamiltonian 
equations of motion and adding some kinetic energy to 
the condensate j = 1 (achieved by the initial condition 
Ppi 7^ 0, which means the excitation of BEG j ~ 1). Fig. 
|4^ shows the time evolution of the radial extension of the 
wave functions in a stack of two BEGs (represented by 




FIG. 5. Normal modes of the coupled radial motion of two 
interacting BECs. The two BECs either oscillate in phase (on 
the left) or with a phase shift of vr (on the right) and different 
oscillation frequencies corresponding to the different normal 
modes. 
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(E)" 



FIG. 6. The frequency Co — — tJ2j/2 that can be inter- 
preted as a characteristic frequency for the energy exchange 
between two condensates of Fig. |4] The frequency depends 
crucially on the distance A between the condensates and the 
scattering length a/ ad and reaches its highest value for small 
distances A and a — )■ arr. 



FIG. 7. Frequency spectra of the radial oscillation of the 
single condensates in the stack of 2 BECs of Fig. [1] for dif- 
ferent excitation energies E*. The single peaks indicate the 
fundamental frequencies, and the peak heights are plotted in 
arbitrary units. 



tigation of these normal modes and the corresponding 
eigenfrequencies shows that the energy exchange signifi- 
cantly depends on the scattering length and the distance 
between the condensates. We demonstrate this for a 2- 
layer stack of BECs: In this case the linearized motion 
can be described by one normal mode where the BECs 
oscillate in phase with a frequency uji and another with 
a phase shift of tt and a frequency uj2 (see Fig. [s]). Thus, 
for initial conditions Upi = 0, Up2 7^ and Upi.2 = the 
deviation from the fixed point is described by 

[ ( + sin(a;2t)\ 
\Up2 J Ysin(wi<) - sin(a;2t) J 

^ /sin(^i)cos(^i)A 
^cos {"^^^t) sin {'^^^t) ) ' 

where Hj = \uii — uj2 \ /2 is the frequency of the envelope 
and can be interpreted as the characteristic frequency of 
the energy exchange. Fig. [6]shows the frequency Cj in de- 
pendence on the scattering length a/ ad and the distance 
A between the condensates. It reaches its highest value 
for small distances and a scattering length near the crit- 
ical value. Increasing the distance as well as increasing 
the scattering length, uj becomes smaller and vanishes 
for A —7> 00 and a /ad 00, respectively. The range of 
frequency shown in Fig. [6] reaches from about w = 1.4 
to a) = 180 which, considering the particle number scal- 
ing, corresponds to about 0.5 Hz to 70 Hz so that this 
effect of energy exchange should be observable in actual 
experiments. 

Our calculations with the exact Hamiltonian equations 
of motion confirm this behavior also for high excitations 
and additionally reveal a dependence of the frequency of 
energy exchange on the excitation energy. In case of high 
excitations the oscillations are anharmonic (see Fig. |4]d) 
and we determine the oscillation frequencies by Fourier 
transforming the time-dependent extension of the BECs. 



Fig. [7] shows the fundamental frequencies of the oscil- 
lations for different excitation energies E*. Additionally 
there appear higher harmonics (not shown) whose am- 
plitudes grow with increasing excitation. The frequency 
of the envelope, however, remains determined solely by 
the difference of the two fundamental frequencies, which 
becomes smaller with increasing excitation energy of the 
stack, i.e. in a highly excited stack the exchange of en- 
ergy between the single BECs takes longer than it does 
in slightly excited stacks. 

We note that, because of the strong confinement in the 
z-direction, the extension of the BECs in this direction 
remains small compared to the distance A also for high 
excitations of the stack so that tunnelling between the 
single condensates can still be neglected. 

This exchange of energy can have drastic consequences 
on the stack of dipolar BECs. We demonstrate this by 
considering two excited BECs. If the condensates are sep- 
arated and their energy is below the saddle point energy 
of the corresponding external potential V the extension 
of the two BECs will oscillate around its stationary value 
for all times. The situation is different if these two con- 
densates are placed in a stack where they interact with 
each other. We show this by exciting two BECs in such 
a way that the energy of each individual BEC is below 
the saddle point energy of V but the excitation energy 
of the whole stack lies slightly above one of the saddle 
points of V. Again, the extensions of both condensates 
oscillate around their stationary values quickly and the 
condensates exchange energy (see Fig. [8^). The last 
few oscillations (see Fig. [8]d) show a highly excited BEC 
(j = 1) while the other one (j = 2) loses its energy and 
settles down in its stationary state. Fig. [HJ; shows the 
BEC with j — 2 that has transferred its whole excitation 
energy onto the first one and resides in its stationary 
state. The whole excitation energy is now located in the 
first condensate whose extension oscillates strongly, until 
at a certain time (t ss 1.0) it becomes so small that the 



8 



a) 0.90 I 1 1 1 1 1 1 r 




0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 

t 




0.960 0.965 0.970 0.975 0.980 0.985 0.990 0.995 1.000 1.005 
t 




0.998 0.999 1.000 1.001 

t 

FIG. 8. Time evolution of the radial extension of the two 
BECs of Fig. |4] described by a particle moving in the external 
potential V with initial momenta p^- ,pzj 7^ and initial values 
for the generalized coordinates that are not the fixed points 
of the Hamiltonian equations of motion. The mean-field en- 
ergy of the system lies slightly above the saddle point of the 
external potential, a) The extension of the BECs oscillates 
and they exchange their excitation energy, b) The last few 
oscillations, c) The second BEG {j — 2) has transferred its 
whole kinetic energy to the first one which is now highly ex- 
cited and whose extension finally reaches qpi — 0, meaning 
the collapse of the condensate. 

attractive contact interaction {a/ad < 0) which is pro- 
portional to the density |7/;(r)|^ becomes dominant and 
the radial extension of the BEC reaches qpi — >■ 0, mean- 
ing its collapse. In the Hamiltonian picture, the particle 
representing the stack of BECs crosses a saddle point 
of the external potential V and falls down at the other 
side. Consequently, excited BECs are able to interact 
with each other in a way that induces the collapse of one 
of the BECs in the stack. Of course this behavior can 
also be observed for more BECs in a stack. 

Finally, we wish to investigate the dynamics of the 
multilayer stack of BECs when the scattering length is 
reduced below the critical value. Above that value, there 
exist several stationary states of the coupled BECs, and 
we assume the stack to be in the stable ground state at 
a scattering length that lies slightly above the critical 
value. If we now decrease the scattering length below 
there no longer exists a stationary state and the resulting 
dynamics of the stack is shown in Fig. [9^ for A = 0.035 
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FIG. 9. Time evolution of a stack of Af, = 3 BEGs. The mul- 
tilayer stack of dipolar BEGs is assumed to be in the ground 
state at a scattering length that lies slightly above the critical 
value and does therefore not change in time (dashed lines). 
Reducing the scattering length below the critical value, all 
coordinates qpj begin to shrink, representing the contraction 
of the BEGs. a) The central BEG (blue line) collapses first 
but at a distance of A = 0.035, this also causes the collapse of 
the other BEGs (dotted line), b) The same situation with a 
distance of A = 0.07. The central condensate again collapses, 
but the coupling is not strong enough to determine the other 
BEGs to collapse. 



for a stack with 3 BECs. The radial "extensions" qpj of 
all BECs begin to shrink until the central BEC {j ~ 2) 
collapses {qp2 = 0). The outer two BECs still have a 
finite extension, but are strongly affected by the collapse 
of the central condensate. The equations of motion are 
solved for the two remaining BECs for the time after the 
collapse of the central BEC (dotted lines) and one can see 
that the other two BECs will also reach qp = 0, meaning 
its collapse. Doubling the distance A between the BECs, 
the interaction becomes weaker. It is again the central 
BEC that collapses first, but the effect on its neighbors 
is not strong enough to force them to collapse, too (see 
Fig. [9]d). 

In the framework of this variational approach, three- 
particle collisions, causing particle losses, have been ne- 
glected, but will inevitably happen during the collapse as 
the density becomes higher. However, we do not expect 
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qualitative changes in the dynamic since the process of 
considering a constant particle number in the BEC dur- 
ing the collapse and afterwards neglecting its influence 
completely (as done here) will then only be changed by 
a continuous decline of the particle number in the cen- 
tral BEC. Nevertheless, we will investigate this by both, 
variational and numerically exact grid calculations, tak- 
ing particle losses into account. 



IV. CONCLUSION AND OUTLOOK 

Describing the multilayer stacks of dipolar BECs varia- 
tionally we were able to show that such a stack is charac- 
terized by several stationary states that arise in tangent 
bifurcations at different values of the scattering length. 
Moreover, we could demonstrate coupled dynamics of the 
stack which reaches from normal modes for small excita- 
tions of the stack and an energy exchange between single 
BECs on experimentally accessible time scales to the in- 



duced collapse of a BEC for high excitations. Generally 
in an excited stack of interacting BECs the individual 
condensates always exchange energy and this exchange 
of energy is significantly affected by the scattering length 
and the distance between the single BECs. 

Of course the ansatz of a single Gaussian trial wave 
function implies restrictions concerning the structure of 
the wave functions of the individual BECs, and even 
though it is not capable of reproducing symmetry break- 
ing angular collapse mechanisms, its advantage has to 
be seen in the description of the dynamics of the stack 
which is easily accessible. To additionally investigate ef- 
fects beyond the Gaussian form of the wave function the 
ansatz can be extended to a Gaussian wave packet for 
each BEC in the stack. In the case of a single BEC, this 
extended ansatz is able to reproduce the numerical re- 
sults [12] what can also be expected for the multilayer 
stack of BECs. Moreover the results found in this paper 
can serve as a useful guide for investigations of the dy- 
namical effects in the stack of interacting BECs by exact 
numerical calculations. 



[1] A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, and 

T. Pfau, Phys. Rev. Lett. 94, 160401 (2005). 
[2] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T. 
Pfau, Rep. Prog. Phys. 72, 126401 (2009). 



[3] S. Ronen, D. C. E. Bortolotti, and J. L. Bohn, |Phys. Rev.| 

Lett., 98, 030406 (2007). 
[4j I. Tik honenkov, B. A. Malomed, and A. Vardi, |Phys.| 

Rev. Lett. 100, 090406 (2008)^ 

"[STO. Dutta and P. Meystre, Phys. Rev. A 75, 053604 

(2007). 

[6] M. Klawunn and L. Santos, Phys. Rev. A 80, 013611 

(2009). 



[7] P. Koberle and G. Wunner, |Phys. Rev. A] 80, 063601 
(2009). 

[8] P. Koberle, H. Cartarius, T. Fabcic, J. Main, and 

G. Wunner, New Journ al of Physics 11, 023017 (2009). 
[9] K. Goral and L. Santos, |PIys. Rev. A 66, 023613 (2002). 
[10] H. Cart arius, T. Fabcic, J. Main, and G. Wunner, PhysT| 
\ Rev. A|78, 013615 (2008). 
[11] A. D. McLachlan, Mol. Phys. 8, 39 (1964). 
[12] S. Rau, J. Main, P. Koberle, and G. Wunner, [Phys. Rev.| 
I A| 81, 031605(R) (2010). 



